/*  scfcasin.cpp  by K.Tsuru  */
// function ID = 9116  since ver 2.18
/*****************************************
SComplex class
It returns arcsin(z).
1) Casin(z) = -i * Clog{i * z + Csqrt(1 - z * z)}.

2) Let z = x [+|-] i*y
arcsin(z) = arcsin(x [+|-] i*y) ([+|-]y >= 0)
= arcsin[ (1/2){sqrt( (1+x)^2 + y^2 ) - sqrt( (1-x)^2 + y^2 )} ]
[-|+]i*arccosh[ (1/2){sqrt( (1+x)^2 + y^2 ) + sqrt( (1-x)^2 + y^2 )} ] (ref. Iwanami II p.224)
= arcsin {(|1+z|-|1-z|)/2} [-|+] i*arccosh {(|1+z|+|1-z|)/2}

When y = 0 i.e. z = x, let arcsin(x) = p, sin(p) = x then
|x| <= 1.0 must be satisfied.
*****************************************/
#ifndef SN_H
#include "sn.h"
#endif

#define CasinZ1Method 1 // Method 2) slightly faster

SComplex Casin(const SComplex& z){
    if(z == 0.0) return SComplex(0.0, 0.0); // arcsin(0) = 0
    if(Im(z) == 0.0){ // pure real z = x
        if( Dabs(Re(z)) > 1.0) z.Real().SetError(z.Real().DOMAIN_ERR, "Casin(z) with z = x(x > 1.0))", 9116);
        return SComplex(Asin(Re(z)), 0.0);  // arcsin(x)
    }
#if CasinZ1Method // 14.1sec
    SComplex zp1(z + 1.0), zm1(z - 1.0);
    SDouble azp1, azm1, rp, ip;

    azp1 = Cabs(zp1);   // |1+z|
    azm1 = Cabs(zm1);   // |1-z|

    rp = DsDiv(azp1 - azm1, 2); // (|1+z|-|1-z|)/2
    rp = Asin(rp);

    ip = DsDiv(azp1 + azm1, 2); // (|1+z|+|1-z|)/2
    ip = Acosh(ip); // ip > 0
    if ( z.Imag().Sign(9115) < 0 ) ip.ChangeSign(9115);

    return SComplex(rp, ip);
#else // 16.4sec
    SComplex p = 1.0 - z * z;
    p = I * z + Csqrt(p);
    p = MI * Clog(p);
    return p;
#endif
}